Skip to content

A new SN feedback model#995

Merged
chongchonghe merged 44 commits intodevelopmentfrom
chong/particles/SN-deposition-new-model
Apr 26, 2025
Merged

A new SN feedback model#995
chongchonghe merged 44 commits intodevelopmentfrom
chong/particles/SN-deposition-new-model

Conversation

@chongchonghe
Copy link
Contributor

@chongchonghe chongchonghe commented Apr 15, 2025

Description

In this PR I implement a new SNe feedback model that smoothly transitions between thermal feedback, kinetic feedback, and momentum feedback, depending on how well the cooling/shell formation phase of the SNR evolution is resolved. An innovation of this scheme is that we do a two stage deposition scheme to account for SNR overlapping and get rid of depedence on the order that SNe events occur.

I have tested three models:

  1. Model 0: Pure thermal deposition (for comparison only, not for production)
  2. Model 1: Thermal energy deposition in the resolved regime, thermal + kinetic in the intermediate regime, and terminal momentum in the unresolved regime.
  3. Model 2: Kinetic energy deposition in the resolved regime, thermal + kinetic in the intermediate regime, and terminal momentum in the unresolved regime.

At a SN event, we distribute mass, momentum, and energy uniformly among the cells that intercect with s sphere of radius 3 dx centered at the center of the cell that the SN is located in. No background smoothing is conducted before deposition.

The results looks great! See figure below. Both model 1 and 2 produces 'correct' and almost identical results. The only difference is the maximum temperature (defined as the maximum temperature in any cell throughout the simulation) in model 2 is slightly lower than that in model 1. In model 2, right after SN explosion, the gas is cool but moves fast, and the shock wave immediately heat the gas to comparable temperature to that in the thermal model. After some experiments, I am convinced that the thermal model and kinetic model (in the resolved regime) has no difference in accuracy or performance related to timestep. This is because, imagine a star explodes into vacuum, q simple calculation will show that the sound speed of the hot SNR in the pure thermal model is about $1/\sqrt{3}$ times the bulk velocity of the cool SNR in the pure kinetic model (microscopic thermal motion has three degrees of freedom while kinetic motion has one). The thermal model turns out to have slightly longer timestep overall. In the TIGRESS model they advocate thermal model because they have over cooling issue, but we don't have this issue at all (thanks to Ben's awesome Grackle cooling module which has a robust ODE solver).

compare

Reproducing Figure 6 of Kim & Ostriker 2017. The spacial resolution and SN ejecta mass/energy/momentum in our simulations are the same as theirs.

Checklist

Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an x inside the square brackets [ ] in the Markdown source below:

  • I have added a description (see above).
  • I have added a link to any related issues (if applicable; see above).
  • I have read the Contributing Guide.
  • I have added tests for any new physics that this PR adds to the code.
  • (For quokka-astro org members) I have manually triggered the GPU tests with the magic comment /azp run.

@dosubot dosubot bot added the size:XL This PR changes 500-999 lines, ignoring generated files. label Apr 15, 2025
@dosubot dosubot bot added enhancement New feature or request particles labels Apr 15, 2025
@chongchonghe
Copy link
Contributor Author

The next step is to test overlapping and non-uniform background. I would need to set up a non-uniform initial condition to mimick realistic ISM. @BenWibking Do you have any idea how to quickly produce such an I.C.?

@BenWibking
Copy link
Collaborator

The next step is to test overlapping and non-uniform background. I would need to set up a non-uniform initial condition to mimick realistic ISM. @BenWibking Do you have any idea how to quickly produce such an I.C.?

I would just run it until you have multiple SN going off in a periodic box, like the RandomBlast problem. Then the background will be very non-uniform :D

Copy link
Contributor

@github-actions github-actions bot left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 10 out of 18. Check the log or trigger a new build to see more.

@chongchonghe chongchonghe changed the title A new SN feedback model [WIP] A new SN feedback model Apr 16, 2025
@BenWibking
Copy link
Collaborator

Regarding model 2, for supernovae that stay confined to the disk, I think it's a fine model. However, I am slightly concerned that this model might artificially inhibit the buoyancy of superbubbles. Buoyancy-driven winds can be significant for some galaxies.

So I think the best test to compare would be of models 1 and 2 in a stratified box, where the Brunt-Vaisala timescale is less than or comparable to the cooling timescale.

@markkrumholz Do you have any thoughts on this particular issue?

@chongchonghe
Copy link
Contributor Author

chongchonghe commented Apr 17, 2025

I've finished a test of our model 2 in the draft, or SNScheme::SN_thermal_or_thermal_momentum in the code, in a more realistic ISM environment. The result is satisfactory: the terminal momentum agrees well with the classical fitting formula at the given density.

The initial condition was created from a uniform background with 27 SN explosions at t=0 and ran till t = 80 Myr to come to equilibrium. SNR overlapping did not cause unexpected behaviour, nor cause the code to crash, thanks to the limiter in our scheme. I'm still testing the accuracy/convergence of our method when SNR overlaps.

CleanShot 2025-04-18 at 02 35 27@2x CleanShot 2025-04-18 at 02 36 01@2x

@BenWibking
Copy link
Collaborator

BenWibking commented Apr 17, 2025

That looks nice!

The only possible improvement that I can think of is to do the deposition using a spherical kernel, such as one of the Wendland kernels (https://link.springer.com/article/10.1007/BF02123482).

There is an implementation of Wendland's $C^2$ kernel already in Quokka here:

AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto kernel_wendland_c2(const amrex::Real r) -> amrex::Real

@chongchonghe
Copy link
Contributor Author

Regarding model 2, for supernovae that stay confined to the disk, I think it's a fine model. However, I am slightly concerned that this model might artificially inhibit the buoyancy of superbubbles. Buoyancy-driven winds can be significant for some galaxies.

So I think the best test to compare would be of models 1 and 2 in a stratified box, where the Brunt-Vaisala timescale is less than or comparable to the cooling timescale.

@markkrumholz Do you have any thoughts on this particular issue?

The model 1 and model 2 in my draft are basically the same model, the only difference being model 2 does a smooth transition from unresolved to resolved phase. Therefore, model 2 is always (slightly) better than model 1. I also observe that in my tests. I don't get how that model could "artificially inhibit the buoyancy of superbubbles".

@BenWibking
Copy link
Collaborator

Regarding model 2, for supernovae that stay confined to the disk, I think it's a fine model. However, I am slightly concerned that this model might artificially inhibit the buoyancy of superbubbles. Buoyancy-driven winds can be significant for some galaxies.
So I think the best test to compare would be of models 1 and 2 in a stratified box, where the Brunt-Vaisala timescale is less than or comparable to the cooling timescale.
@markkrumholz Do you have any thoughts on this particular issue?

The model 1 and model 2 in my draft are basically the same model, the only difference being model 2 does a smooth transition from unresolved to resolved phase. Therefore, model 2 is always (slightly) better than model 1. I also observe that in my tests. I don't get how that model could "artificially inhibit the buoyancy of superbubbles".

Your description of the models was:

Model 1: Thermal energy deposition in the resolved regime, thermal + kinetic in the intermediate regime, and terminal momentum in the unresolved regime.
Model 2: Kinetic energy deposition in the resolved regime, thermal + kinetic in the intermediate regime, and terminal momentum in the unresolved regime.

So if you only inject kinetic energy in the resolved regime, then you are not directly increasing the entropy of the gas. The entropy of the gas will increase after it shocks, but it may be delayed compared to when the ISM entropy should increase, so any high-entropy bubbles may not rise as far as they should.

@BenWibking
Copy link
Collaborator

Your plot of the maximum temperature with this model suggests there could be an effect like this, since the maximum temperature is less for Model 2 compared to Model 1 in the resolved phase.

@chongchonghe
Copy link
Contributor Author

Your plot of the maximum temperature with this model suggests there could be an effect like this, since the maximum temperature is less for Model 2 compared to Model 1 in the resolved phase.

Oh, I was referring to a different 'model 1' and 'model 2'. Never mind. I will give a thorough reply tomorrow.

@chongchonghe
Copy link
Contributor Author

chongchonghe commented Apr 18, 2025

I am now convinced that model 1, "Thermal energy deposition in the resolved regime, ...", is as good as model 2, as long as we have a reliable cooling module like the one in Quokka. There is no difference in timesteps between the two models, as I argued above. The only advantage of model 2 occurs when SNe are ejected into vacuum and the temperature goes to 10^9 K; an inaccurate cooling module (e.g. a single Backward Euler) could lead to overcooling or an inaccurate cooling curve. Besides, while both models have similar timesteps, model 1 places more pressure on the cooling step, which could negatively impact performance if it requires many iterations.

model 1 and model2 are only different when RM < 0.027, which occurs at n <~ 0.1 cm-3. The realistic ISM test I did above has a background density of ~0.5 cm-3, so these two models make no difference. Considering this, I reckon model 2 does nothing to prevent superbubbles, because pure kinetic deposition is only activated when the shell formation is well resolved (dx < 0.1 R_sf. Note that we use 3 cells and 0.3^3 = 0.027), and in this regime, the ambient gas gets immediately shock-heated and the subsequent evolution is very similar between the two models. See temporal profiles below and pay attention to the last panel -- the max temperature v.s. time in log-log scale.


Fig1: Temporal profile of the Sedov-Taylor phase of a SN event at a background density of n = 0.1, using model 1. Fig2: same as Fig1 but with model 2.

@chongchonghe
Copy link
Contributor Author

/azp run

@azure-pipelines
Copy link

Azure Pipelines successfully started running 2 pipeline(s).

@chongchonghe
Copy link
Contributor Author

@BenWibking I tried on my computer and all the tests that use GrackleLikeCooling failed with SIGILL Invalid caught by amrex.fpe_trap_invalid=1. Can you take a look?

I think we really need to add one of the test problems with GrackleLikeCooling into the unit tests. Currently none of those problems are activated. We can reduce max_timesteps to make ti run fast. That will help catch any problem caused by Grackle cooling module.

@chongchonghe
Copy link
Contributor Author

chongchonghe commented Apr 22, 2025

I have to temporarily disable FPE in test_SN.

Another option is to disable Grackle cooling for test_SN.cpp. The SN feedback module does not rely on cooling anyway.

@chongchonghe
Copy link
Contributor Author

/azp run

@azure-pipelines
Copy link

Azure Pipelines successfully started running 2 pipeline(s).

github-merge-queue bot pushed a commit that referenced this pull request Apr 23, 2025
### Description

The previous `Gravity3D` problem tested too many things simultaneously,
making it hard to track specific module functions. Consequently, I have
replaced it with a new test called `ParticleCreation`, which tests
particle creation from cell, stellar evolution, and particle
destruction. It's unnecessary to test gravity in this problem because
gravity is tested in `BinaryOrbit`. Similarly, SN feedback is tested in
the `SN` test in #995. Additionally, `Gravity3D` is badly named.

More changes:

- Add logging of the number of particles created inside
`createParticlesImpl`
- Add logging of the number of particles removed at each level inside
`destroyParticlesImpl`
- `updateEvolutionStage` updates a particle's stage based on its
`deathtime`.
- Add an extra evolution stage: `Removed`. Only those particles will be
removed.
- Updated the format of `printParticleStatistics` 

### Checklist
_Before this pull request can be reviewed, all of these tasks should be
completed. Denote completed tasks with an `x` inside the square brackets
`[ ]` in the Markdown source below:_
- [x] I have added a description (see above).
- [x] I have added a link to any related issues (if applicable; see
above).
- [x] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [x] I have added tests for any new physics that this PR adds to the
code.
- [x] *(For quokka-astro org members)* I have manually triggered the GPU
tests with the magic comment `/azp run`.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
@dosubot dosubot bot added size:XL This PR changes 500-999 lines, ignoring generated files. and removed size:XXL This PR changes 1000+ lines, ignoring generated files. labels Apr 24, 2025
@sonarqubecloud
Copy link

@chongchonghe
Copy link
Contributor Author

/azp run

@azure-pipelines
Copy link

Azure Pipelines successfully started running 2 pipeline(s).

@chongchonghe
Copy link
Contributor Author

This PR is well tested and should work just fine. @BenWibking If you want to use this for a test run for your talk or proposal, go for it.

@BenWibking
Copy link
Collaborator

This PR is well tested and should work just fine. @BenWibking If you want to use this for a test run for your talk or proposal, go for it.

A star formation model to create the star particles is still needed though, right?

@chongchonghe
Copy link
Contributor Author

This PR is well tested and should work just fine. @BenWibking If you want to use this for a test run for your talk or proposal, go for it.

A star formation model to create the star particles is still needed though, right?

Right. We'll need Aditi's StochasticStellarPop module.

@markkrumholz
Copy link
Collaborator

Just want to check on the status of this code, since I'm now back online. Is there stuff I should review, or @BenWibking are you happy to do so? I will also go look at the paper draft in any event.

@chongchonghe
Copy link
Contributor Author

chongchonghe commented Apr 26, 2025

Just want to check on the status of this code, since I'm now back online. Is there stuff I should review, or @BenWibking are you happy to do so? I will also go look at the paper draft in any event.

@markkrumholz Oh. I specifically requested waiting until your review before we merge this branch. Ben has reviewed it and I have addressed all his comments.

Copy link
Collaborator

@BenWibking BenWibking left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me.

@dosubot dosubot bot added the lgtm This PR has been approved by a maintainer label Apr 26, 2025
@chongchonghe chongchonghe added this pull request to the merge queue Apr 26, 2025
Merged via the queue into development with commit 50c179c Apr 26, 2025
24 of 25 checks passed
@chongchonghe chongchonghe deleted the chong/particles/SN-deposition-new-model branch May 12, 2025 06:22
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request lgtm This PR has been approved by a maintainer particles size:XL This PR changes 500-999 lines, ignoring generated files.

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants